function [H, H1] = gauss_fker(siz, a, b, angle)
%GAUSS_FKER Gaussian Filtering Kernel
%
%   A Gaussian kernel is defined to be 
%
%       h(x, y) = exp( - (x - x0)^2 / (2 * a^2) - (y - y0)^2 / (2 * b^2) )
%   
%
%   H = GAUSS_FKER(siz, a);
%   H = GAUSS_FKER(siz, a, b);
%
%       Constructs a Gaussian filtering kernel H.
%
%       Input arguments:
%       - siz:      The kernel siz, can be scalar (if it is square)
%                   or in form of [h w].
%
%       - a:        The standard deviation along the x-axis
%
%       - b:        The standard deviation along the y-axis. 
%                   (If b is omitted, it assumes a = b).
%
%
%   H = GAUSS_FKER(siz, a, b, angle);
%
%       Constructs a rotated Gaussian kernel. In particular, it returns
%       the rotated version of GAUSS_FKER(siz, a, b) by the specified
%       angle (counterclockwisely).
%
%   [H, H1] = GAUSS_FKER(siz, a, b, angle);
%
%       Additionally constructs the (1st order) Derivative of Gaussian 
%       filter, along the specified direction.
%
%

% Created by Dahua Lin, on April 3, 2012
%

%% parse arguments

if isscalar(siz)
    h = siz;
    w = siz;
elseif numel(siz) == 2
    h = siz(1);
    w = siz(2);
else
    error('gauss_fker:invalidarg', 'siz should be either a scalar or a pair.');
end

if ~(rem(h, 2) == 1 && rem(w, 2) == 1)
    error('gauss_fker:invalidarg', 'kernel dimensions should be odd integers.');
end

if nargin < 3
    b = a;
end

if nargin < 4
    angle = 0;
end

%% main

[X, Y] = fker_xymap(h, w, angle);
H = exp( - (X.^2) * (1/(2*a^2)) - (Y.^2) * (1/(2*b^2)) );

if nargout >= 2
    H1 = H .* (X * (-1/a^2));
end



